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ABSTRACT 



The magnetorotational instability (MRI) is believed to be responsible for most of the angular momentum transport in accretion discs. 
However, molecular dissipation processes may drastically change the efficiency of MRI turbulence in realistic astrophysical situations. 
The physical origin of this dependency is still poorly understood as linear and quasi linear theories fail to explain it. In this paper, 
we look for the link between molecular dissipation processes and MRI transport of angular momentum in non stratified shearing box 
simulations including a mean vertical field. We show that magnetic helicity is unimportant in the model we consider. We perform 
a spectral analysis on the simulations tracking energy exchanges in spectral space when turbulence is fully developed. We find that 
the energy exchanges are essentially direct (from large to small scale) whereas some non linear interactions appear to be non local 
in spectral space. We speculate that these non local interactions are responsible for the correlation between turbulent transport and 
molecular dissipation. We argue that this correlation should then disappear when a significant scale separation is achieved and we 
discuss several methods by which one can test this hypothesis. 
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1. Introduction 

The transport of angular momentum in astrophysical discs is 
a central problem of accretion theory. To explain discs life- 
time and accretion rates, it is often assumed that these objects 
are turbulent. Turbulence is then included in global models us- 
ing a turbulent viscosity p rescription as in the a disc model 
dShakura & Sunvaevll 1973b . 

The origin of this turbulence has been the subject of 
many debates over the past decades. It is now generally as- 
sumed that the magnetorotational instability, or shortl y MR I 
dWukhov 1959; Chandrasekhar 1960; Balbus & Hawlevi ll99lh . 
is responsible for disc turbulence, although hydrod y namic 
processes might also be at work dLesur & Qgilvid 120 lot 
iLesur & Papaloizoul l2010h . Although MRI generated turbu- 
lence is gene rally efficient at transporting angular momentum 
(iHawlev et al.i,1995) . recent results have shown a strong sensi- 
tivity of MRI turbulence on small scale dissip ation processes 
dLesur & Longaret"till2007l; iFromang et al.ll2007 ). and in partic- 
ular on the magnetic Prandtl number Pm (ratio of microscopic 
viscosity to resistivity). This effect, called the a-Pm correlation, 
casts doubts on the actual efficiency of the MRI in realistic situ- 
ations since Pm can va ry by several order of magnitude in discs 
dBalbus & Henri"2008). Several attempts have been made to ex- 
plain this correlation, either from the linear theory of dissipative 
MRI modes (Pessah & Chan 2008) or fr om the quasi-linear par- 
asitic modes theory (.Pessah & Goodnianll2009l) . However, these 
approaches were shown to be unsu ccessful when compared t o 
high Reynolds number simulations (Longarett i & Lesunl2010h . 
Instead, Longaretti & Lesur (2010) suggested that the a - Pm 
correlation could be due to the nature of the MHD cascade in 
MRI generated turbulence, in which one might expect inverse 



cascades and/or non local interaction in spectral space. This 
kind of process would allow for a direct communication between 
the injection scales (transport scales) and the largest dissipation 
scale (either resistive or viscous). 

The purpose of the pres ent work is to inves t igate some of 
the conjectures presented bv iLongaretti & Lesuil d2010h regard- 
ing the nature of the MHD turbulent cascade in accretion discs. 
To this end, we co nsider several of the high resolution simu- 
lations presented by ILongaretti & Lesti^ d201Q) and we analyse 
the energy exchanges in spectral space. This paper is organized 
as follows. We describe our model, equations and the spectral 
analysis we use in section |2l Section [3] is the core of this paper 
and discusses our numerical results. The implications of these 
results are presented in the final section. 

2. Model and spectral analysis 

2.1. Equations 

In the following, we will adopt the shearing box model which 
accu rately represent the local physics of an acc retion disc 
(Ha wlev et all 119951: [Balbus 2003; Regev^&Umurhan 20081. 
The adopted coordinate system is such that x - (r - ro) and 
y = ro(p where ro is the fiducial radius of the shearing box in the 
disc, and (p the azimuthal coordinate in the rotating frame. The 
velocity can be decomposed as a mean velocity plus a fluctuat- 
ing part U = -qQ.xey + v where Q is the local rotation rate and 
q - 3/2 for a Keplerian disc rotation profile. As a simplifica- 
tion, we assume the flow is incomp ressible, which correspond s 
to the "small shearing box limit" of lUmurhan & RegevI d2004l) . 
The equations of motion then read 

d,v = -V - Vv -Vn + B -VB + q^lxdyv 
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+2Q.Vyex - (2 - q)nv_^-ey + vV v, 
d,B = -V VB + B Vv 

+qQ.xdyB - qQ.B,ey + rjV^B, 
V • v = 0, 
V • B = 0, 



(1) 
(2) 



where TL is the total pressure, v the viscosity and r/ the ohmic 
resistivity. In the following, we impose a mean vertical field Bq 
which will be conserved during the evolution of the flow thanks 
to the shearing-sheet boundary conditions. It should be noted 
that the magnetic field strength is expressed in Alfven speed for 
simplicity. 

Several dimensionless numbers characterise the equations of 
motions. In this paper, we will use the following: 

- The amplitude of the imposed mean vertical field measured 
by 



y8 = 



(3) 



where L, is the vertical box size. This definition mimics the 
usual plasma (3 in vertically stratified discs obeying the ver- 
tical hydrostatic equilibrium constraint Cj ~ Q.L^ 
- The viscous Elsasser number: 



Av 



K 

Qv 



(4) 



which is related to th e Reynolds number used in 
iLongaretti & Lesud (l2010l) by A^ = qRelp 
- The resistive Elsasser number: 



A„ 



Q.r] 



(5) 



connected to the magnetic Reynolds number by a similar re- 
lation 

- The magnetic Prandtl number 



Pm 



A, 
Av 



(6) 



which compares the amount of viscous and resistive dissipa- 
tion. 



2.2. Fourier transform in slieared flows 

It is convenient to introduce the shearing frame (x', y', z') : 

X — x' , 

y - y' - qCltx', 
z - z ■ 

Writing the equations of motions in the sheared frame allows us 
to eliminate the explicit spatial dependency: 



d,v = -v-Vv- Vn + B ■ VB 

+2Q.Vyex - (2 - q)Q.Vxey + vV^v, 

d,B = -V VB + B Vv 
-qQ.B_,ey + rjV^B. 



(7) 



(8) 



where v and B are now assumed to be functions of x' so that 
the nabla operator expression becomes 



In the sheared frame, the shearing sheet boundary conditions 
make every physical quantity X(x') periodic so that X can be 
expanded in Fourier series: 

X{x') = ^Xfc' expOfc' ■x')^Yj ^Ht) exp(/fc(0 ■ x) . (10) 

k' k(t) 

This relation defines the time dependent unsheared wave vectors: 



kx - kx' + qQ.ky't, 

ky 

k. 



k,,. 



(11) 
(12) 
(13) 



As expected, the application of the V operator to X{x') corre- 
sponds to a multiplication of its Fourier components by ik{t). 

This definition of the unsheared wave vectors with the 
Fourier decomposition ( fTOl ) is usually refeiTed to as a "shear- 
ing wave" decomposition of a sheared flow. It was first 
used by Lord K elvin to study the stability of sheared flows 
(Thomson '1887) and later in the astrophysical context by 
Goldreich & Lvnden-Bell (.1965.) for spiral arms of galaxies. 

2.3. Sliell filter decomposition 

Following lFrischI (Il995h and lAlexakis et al.l (l2007b . one defines 
shell filtered quantities in the unsheared Fourier space. At any 
given time, a series of linearly spaced shell sizes K is defined 
from Ki to ^T^ax; by construction 6K - Kj - Kj^x is the shell 
width (for any j). We define the shell-filtered field Xk in shell 
Kjhy. 



Xk 



^ Xk exp (ik 

Kj-SK/2<k<Kj+SK/2 



X). 



(14) 



(9) 



It should be noted that since k depends on time, the exact num- 
ber and the distribution of the modes entering the above formula 
for any given K might vary in time, adding an extra complicatio n 
compared to the homogeneous case of Alexakis et al ] (l2Q07h . 
This point is discussed in the Appendices. 

2.4. Energy transfer equations 

The transfers we are interested in relate to the equations of 
the kinetic and magnetic energies. These shell-restricted, box- 
averaged equations involve a number of transfer functions that 
are introduced along with the related equations. We follow here 
the logic of Alexakis et al. (2007) and extend it to shear flows. 
The analysis of transfers gives indications about the locality of 
interactions in Fourier space. 

Because of the incompressibility condition, energy transfers 
(but also stress and magnetic helicity) involve only at most the 
product of three components of the velocity and magnetic fields 
and their derivatives. Therefore, couplings in Fourier space de- 
pend only on triads of wave-vectors, noted fci, fc2 ks- The clos- 
ing condition (fci + k2 - ks) furthermore imposes that at least 
two of the wave vectors ki are of the same magnitude; the third 
one can be either much smaller (implying nonlocal couplings 
through large scales) or of the same magnitude as the other two 
(implying locality of couplings in Fourier space). 

Using the equations of motion in the sheared frame (|7]— [8]l, 
one can derive the equation for the shell filtered energy density: 

d,EK = 2 [^"'^2, K) + TUQ, K)] + S ,,K 
Q 
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+qQ.Iy,K - vDy^K 

d,MK = ^ [Tbb(Q,K) + T,b(Q,K)] + S h.x 

Q 

+qQ.Ib,K - rjDh^K 



(15) 



(16) 



where Eg = and Mk - {B^/2}. In the above expression, 

we have defined the transfer functions by 



TUQ,K) = -{vk-(v-V)vq), 
Tbb(Q,K) = -{Bk-{v-V)Bq), 
Tbv(Q,K) = {vk-(B-V)Bq), 
T.b(Q,K) = {Bk-(B-V)vq), 

where ( . ) denotes an spatial average on the shearing box vol- 
ume. Tij(Q,K) represents the transfer from energy (kinetic 
or magnetic) from shell Q to energy '/ (kinetic or magnetic) 
in shell K. Note that Ti/Q, K) - -Tji(K, Q), so that whatever is 
taken from one shell of one type is totally transfeiTed to the other 
shell. Similarly, TaiK, K) - 0: there is no effective transfer from 
a shell to itself, as should be. This justifies the identification of 
these quantities as shell-to-shell energy transfer functions; in ef- 
fect, the third member of a triad is only a relay in eff'ective energy 
exchang es between sh ells Q and K (for a more detailed discus- 
sion, see IVermall2004 . In these expressions, Vk ■ Vk - Vk ■ v, 
has been used, as well as v - Yjq vq. 

The next terms in the shell energy budget involve energy 
transfers due to the mean shear Sr,K and S h,K- These terms are 
singular in time as they correspond to energy fluctuations due to 
wave entering or leaving the shell K as k{t) evolves. They can be 
formally defined by: 



Sx.K = 2_i — 2 — ~ *' 



= z 



qQ.kyk_^{t) X*j^,Xk' 

m\ ~ 



x\^6{\k{t)\ -K + 6KI2^ - 5^ - \k(t)\ + K + dKIl^ 



where tk' is the instant when the wave fe' enters or exits the shell 
K and e^- = +1 for an entering/exiting wave (see appendix lAllBl ). 

The remaining terms in the shell kinetic and magnetic energy 
budgets, 

Iv,K = {v.x,KVy,K), 
h,K - -{Bx,KBy^K), 
D,,K = <(V X Vk)\ 

Dh,K = <(V X Bk?}, 

represent the energy injection through the shear and dissipation 
through viscosity and resistivity in shell K. 

2.5. Energy fluxes in spectral space 

Using the transfer function defined above, it is possible to intro- 
duce energy flux in Fourier space. In this work, we will use the 
following fluxes: 



rh(Ko,t) = |] 2rw,(e,/:) 

K=Ko Q 

TAKo, 0=22 K) + TUQ, K) 



K=Ko Q 



- ^^(1^(01-^) 

ZqQk,,kAi) B*,, ■ B^' , \ 



(18) 

(19) 
(20) 
(21) 



where IF,,, Fi, and 'Fj are respectively the kinetic, magnetic, 
exchange and shear fluxes; the shear fluxes are evaluated in 
K - Ko + dK/2. Each of these fluxes computes the amount the 
energy transferred from shells K < Ko to shells K > Ko (i.e. the 
flux of energy "through" the outer boundary of shell Ko). The ki- 
netic (respectively magnetic) fluxes compute transfers of kinetic 
to kinetic (respectively magnetic to magnetic) energy, whereas 
the exchange flux is a flux of total energy (magnetic plus ki- 
netic) in which magnetic and kinetic energy are constantly trans- 
formed into one another. Shear fluxes are singular in time as they 
are non zero only when a wave enters or leave shells K > Kq. 
It should be noted that shear fluxes are statistically non zero 
only for anisotropic turbulence, as the amplitude of any mode 
(kjc, ky, k.) should be statistically equal to the amplitude of the 
mode (-kx, ky, k^) in isotropic turbulence. 

These fluxes allow one to check the direction of the en- 
ergy transfers due to the nonlinear terms. Indeed, a direct en- 
ergy transfer (large to small scale) implies a positive flux with 
the above definition, whereas an inverse cascade can be charac- 
terised by a negative flux. 

2.6. Numerical method 

Equation ([1]) and (|2| are solved using the Snoopy code. Snoopy 
is a 3D spectral (Fourier) method based on the shearing wave 
decomposition (fTTT l — (fT3] l. The Fourier transforms are evalu- 
ated using the FFTW 3 Ubrary, with both MPI and OpenMP 
parallelisation techniques. N onlinear terms are c omputed using 
a pseudo-spectral algorithm (ICanuto et al.lll988h and antialias- 
ing is enforced using the "3/2" rule. Time integration is per- 
formed by a third order, low storage Runge-Kutta scheme for 
non-linear terms, whereas an implicit scheme is used for vis- 
cous and resistive terms. This spectral scheme uses a periodic 
remap algorithm in order to continuously follow the smallest 
wav e number of the system in the sheared frame |fc(f)l < ^max 
(see lUmurhan & RegevI |2004| appendix C for a complete de- 
scription of the periodic remap algorithm). Note that the periodic 
remap method used in this code i s diffe rent from the continuous 
remap method used bv lLithwickldloOTi) . Our main motivation to 
implement a periodic remap is the possibility of using the 3/2 an- 
tialiasing rule and power of 2 grid sizeflfor which Fourier trans- 
form and parallelisation methods are more efficient. This code or 
its variant has been used i n several context, including the MRI 
(iLongaretti & Lesu3l2010h and the subcritical baroclinic insta- 



(17) 



K=Ko Q 



' If the number of point in one direction is a multiple of 2, waves at 
the Nyquist frequency do not have any properly defined phase. This is 
not a problem for classical spectral methods or for the periodic remap 
method since the Nyquist frequency is either in the dissipation range or 
in the antialiasing dump zone. However, when using a continuous remap 
method, the Nyquist frequency waves are remapped to large scale waves 
in physical space, which might lead to non-physical behaviours. 
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Fig. 1. Energy spectrum at Pm - 0.0625 (left) and Pm - 0.25 (right). In the Pm = 0.25 case, a power law spectrum is observed for 
the kinetic energy corresponding to a k^^^^ spectrum. 
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Fig. 2. Energy Injection spectrum at Pm = 0.0625 (left) and Pm = 0.25 (right). Although the injection is significantly reduced at 
small Pm, shape of the spectrum is similar and dominated by the largest scale. 



biUty dLesur & Papaloizoull2010l) . It is available for download 
on the author's website. 



3. Results 

3.1. Simulations parameters and averaging procedure 

The spectra and transfers presented in this section are all derived 
from two simulations of an MRI saturated state. These runs cor- 
responds to the P m = 0.25 and Pm = 0. 0625 high resolution 
runs discussed in iLongarett i & Lesur (2010). Both runs have a 
resolutioiflA^.i xA^^ xA^2 = 768 x 384 x 192 with a box aspect ra- 
tio LxXLyXL^ =4x4x1. We impose a mean vertical field in the 
box withy8 = 10^ and A„ = 30 (/?e = 2 x 10"*). Each simulation is 
integrated for 50 orbits starting from random noise and the spec- 
tra are averaged from the last 40 orbits to remove any influence 
of the initial conditions. The two simulations considered in this 
section only differ by their ohmic resistivity, the Pm = 0.25 run 
having A^, = 7.5 {Rm - 5000) whereas the Pm = 0.0625 run has 
A, = 1.87 (Rm = 1250). 



^ The quoted resolution now includes the aliasing domain. This con- 
vention differs from the one adopted in our previous papers, where only 
the "useful" domain was accounted for when quoting resolutions. 



The statistical average of any quantity X of interest {X)stat 
should in principle be computed on different realizations but are 
evaluated in practice as usual via an ergodic hypothesis: 



X{t)dt^^{ X(t)dt^yX(Td, (22) 

7^^°° ^ Jo ^ Jo ; 



where < Tj < T are a sufficiently large number of instants of 
flow snapshots. 

The spectra are averaged in the spherical shells introduced 
in (fT4l i. The shells are defined so that K„ = 2nn/L, and 6K = 
2nlL.. This means that some power is present in the shell K -Q, 
as it contains large scale horizontal waves with no vertical struc- 
ture. The shell-integrated spectra and transfers obtained by this 
procedure are then averaged in time over 40 instantaneous snap- 
shots (1 snapshot per orbit). For simplicity, we have renormal- 
ized the wave vectors K so that K' - K/2n on all the plots in 
this section. Note that shells K > 32 are incomplete in the y 
direction since the resolution per scale height is smaller in that 
direction. This is not a problem since these shells are in the dis- 
sipative range and high ky modes are weaker than the equivalent 
high k jk~ modes due to the anisotropy of MRI turbulence (see 
section l3T2l i. Note that, with the procedure we have used, one can 
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reconstruct the box averaged quantities by summing the spectra 
over the integers K. 

The shear transfer terms and Fs,h are computed in a spe- 
cial way. Indeed, one cannot compute Fs for a given shell and 
snapshot time numerically because of the 6 functions. Instead, 
we introduce a shell-averaged flux: 

, r^Ko+6KI2 

F-^(Ko,t)=— dKTsiKj). (23) 

"'^ JKo-SK/2 

As {Ts)stat depend only on K (the turbulence is statistically 
stationary) and varies little with K on scales of the order of 6K, 
one has {'T's}siai -{F';^}srai- One can therefore use in the aver- 
aging procedure described above to estimate l^^stat- 

The numerical flux we obtain is then averaged over time fol- 
lowing the procedure described above. 

3.2. Spectra and energy injection 

We first present the energy spectra on Fig. [T] for Pm = 0.0625 
and Pm - 0.25. The standard deviation, measured from 40 in- 
stantaneous snapshots, is shown as a shaded region on these 
spectra. This dispersion is due to temporal fluctuations of the 
turbulence intensity. The most obvious feature observed in these 
spectra is the presence of a k'^^~ spectrum for the kinetic energy; 
the traditional Kolmogorov scaling k^^^^ appears to be excluded 
in the Pm = 0.25 run, but it cannot be strictly excluded in the 
Pm - 0.0625 run. A k^^^^ power-law was also found in zero 
net flux MRI turbulence (Fromang 2010), although the spec- 
trum shape differs, our spectra being exempt of any "bump" at 
intermediate scale. As our runs do not yet resolve the inertial 
range of the turbulent cascade, these apparent spectral shapes 
require some comment. The presence of ajc^ spectrum is usu- 
ally related to the theoretical argument of^roshnikov (1963) and 
iKraichnali (i 19651) (or shortly IK). However MRI turbulence is 
not strongly magnetized, and therefore falls outside the domain 
of validity of the IK phenomenology. Moreover, the magnetic 
field spectrum does not follow any well-defined power law, as 
expected from the wide and overlapping injection (see below) 
and dissipation spectra, indicating that the spectrum we get is 
not an inertial spectrum. We are therefore forced to conclude that 
although the kinetic spectrum looks like an IK or Kolmogorov 
spectrum, it is not described by the I K or K olmogorov phe- 
nomenologies, nor by recent extensions iBoldvrev 2005). 

Changing the magnetic Prandtl number does not change the 
power-law index for the kinetic energy. We note however 2 major 
effects: the overall spectra amplitudes are reduced and the dissi- 
pation scales move to larger scale as one reduces Pm. These two 
effects are expected since it is known that smaller Pm turbulence 
is associated with a smaller transport efficiency and therefore a 
weaker injection of energy in the cascade. This effect is con- 
firmed by the injection spectra (Fig. |2]i which are significantly 
reduced at smaller Pm. 

We note that the energy injection peaks at the largest scale 
of the box, although injection still exist at A: ~ 10. Therefore, 
although a power-law spectrum is found for 2 < A; < 10, this 
spectrum cannot be described as an "inertial range" since energy 
is still injected at these intermediate scales. 

We present on Fig.[3]bidimensional spectra of magnetic en- 
ergy for Pm = 0.25 (kinetic spectra are not shown as they share 
essentially the same properties). These spectra were obtained av- 
eraging 3D energy spectra over 40 orbits and taking the average 
in the k^, ky, or k, directions. We first note a strong anisotropy 
in the (kjc, k^) plane which indicates that trailing shearing waves 



{kxky > 0) have more energy than leading shearing waves 
{k^ky < 0). As we will see below, this results in non zero shear 
transfer terms. 

Looking at the aspect ratio of the energy contours, we see 
that turbulence is slightly less anisotropic at large k than at small 
k (the contours are less "elongated" at large k), although com- 
plete isotropy is not yet reached in this simulation. Let us how- 
ever remark that the spectral truncation (due to the finite res- 
olution) tends to deform the contours at large k, which might 
accelerate the return to isotropy. One should therefore perform 
higher resolution runs (or at least double Ny) in order to confirm 
this return to isotropy. In principle, one would expect a return 
to isotropy at small scales if the nonlinear transfer terms domi- 
nate all the other terms (injection, body forces) at large enough 
k. However this is not always the case (e.g. in the presence of a 
strong mean magnetic field). 

The {kji,k,) spectrum shows that turbulence is essentially 
isotropic at large k in that plane. For k ~ 1, we find a slight 
anisotropy where modes with k. + are favoured. This is prob- 
ably a result of large scale MRI unstable modes which all have 
k~ in the presence of a mean vertical field. Note that this 
anisotropy disappears very quickly as one moves to larger k. 

3.3. Magnetic lielicity and cross lielicity 

The presence of kinetic and/or magnetic helicity in MHD tur- 
bulence is often invoked to explain large scale dynamo ac- 
tion. Indeed, it is known that an inverse cascade of mag- 
netic heli city can appear in fully developed helical MHD tur- 
bulence (iFrisch et al.l Il975l) . potentially leading to a buildup 
of magnetic energy at large scale. Some kinematic ef fects of- 
ten used in mean field dynamos, like the a effect (iMoflattl 
il978D . also lead to the g eneration to large scale helical 
fields (iBrandenburg & Subra manian 2005). Magnetic helicity 
has therefore been suggested as a possible driving mecha - 
nism (or at least a tracer) of disc dynamos (iBlackmanI |20 1 Oh . 
Moreover, several authors have tried to link magnetic helicity 
conservation and magnetic helici ty flux to t he saturation proper- 
ties of the MRI CVishniac 2009; Kapyla & Korpi 2010). 

In this work, we define the magnetic helicity "H^ = {A-{B- 
(-B))), where A is the vector potential of the fluctuation (this 
expression is gauge-invariant in the shearing box); (•> denotes a 
volume average. We show on Fig.|4](left) the spectrum of relative 
helicity K\Hk\I2Mk for the Pm - 0.25 run. As it can be seen, the 
relative helicity is less than 1% for all scales of these simulations. 
Moreover, this quantity is strongly fluctuating and its sign is not 
well definecfl 

These results tend to indicate that magnetic helicity is dy- 
namically unimportant in the unstratified simulations presented 
here and that MRI saturation is not related to a magnetic quench- 
ing effect due to magnetic helicity accumulated at large scale. 
This was to be expected in the first place as unstratified shearing 
boxes (both with and without mean field) are mirror symmetric. 
Note however that this picture might change when stratification 
is included. 

An other quantity of interest which might play a role in the 
MHD tu rbulence cascade is the cross helicity = {v • b) 
(see e.g. iPerez & Boldvrevll2010l) . When cross helicity is non 
zero, the energy of Alfven waves traveling along and against 
the mean field are not equal. For this reason, turbulence with 
cross helicity is often called imbalanced turbulence. Locally im- 

^ Note that the absolute value of the relative magnetic helicity is plot- 
ted on Fig.|4] 
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Fig. 3. Bidimensional magnetic energy spectra at Pm = 0.25. Left: (k^;k^) spectrum, averaged in ky. Centre: {kx,ky) spectrum, 
averaged in k,. Right: (ky, k,) spectrum averaged in kx- Each contour correspond to a factor 10 in magnetic energy. 





Fig. 4. Average relative magnetic heUcity (left) and cross helicity (right) spectra in the Pm - 0.25 case (black lines). Instantaneous 
spectra are represented in light blue. Note that the absolute value of relative helicities is plotted here, as helicities sign is constantly 
changing. 
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Fig. 5. Energy fluxes as computed from equation (fT7ll-(l2ni at Pm - 0.0625 (left) and Pm - 0.25 (right). Energy fluxes are always 
direct (from large to small scale) and dominated by magnetic and exchange fluxes. 



balanced turbulence is often observed in strong MHD turbu- ulations [Fig. |4](right)]. As for the magnetic helicity, we find that 

lence, where the guide field is weaker than the turbulent flue- the relative cross helicity is small (< 10 ') and highly fluctuat- 

tuations. To check whether non stratified MRI turbulence was ing at all scales. This tends to indicate that MRI turbulence is 

imbalanced, we have computed cross helicity spectra of our sim- not imbalanced in our setup. As for magnetic helicity, this result 
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was to be expected because of the mirror symmetry properties 
of the non stratified shearing box. The absence of any significant 
cross heHcity also shows that energy spectra in Elsasser variables 
= V + are equal and proportional to the kinetic plus magnetic 
energy spectrum. 

3.4. Energy fluxes 

The energy fluxes (fT7])-(l2T]) allow one to check the average di- 
rection of the energy flux in spectral space. To explain the de- 
pendence of the turbulent transport on Pm, seve ral authors (e.g. 
iLesur & Longarettill2007l; iFromang et al.ll2007l) have suggested 
that an inverse cascade driven by resistive and viscous scales 
might be at work. Since magnetic helicity is irrelevant to this 
problem, only the kinetic, magnetic, exchange and shear fluxes 
are important for the non stratified shearing box and should be 
checked for an inverse cascade. 

We present the energy fluxes at Pm - 0.25 and Pm - 0.0625 
on Fig. |5] Standard deviations are shown for kinetic and mag- 
netic fluxes as shaded regions. These deviations are computed 
following the procedure described in section [J!2l We always find 
positive fluxes, meaning that the non linear transfers sae forward 
or direct (from large to small scales) on average. However, at 
larger scale, the standard deviation may allow for an inverse cas- 
cade of kinetic energy in the Pm = 0.25 run. This indicates 
that, although the kinetic cascade is direct on (time) average 
over most of the spectrum, inverse cascades can sometimes be 
observed on the largest scales. This inverse cascade of kinetic 
energy could be an explanation to the large scale hydrodynamic 
structures which are obse rved in several MRI turbu lence sim- 
ulatio ns, such as vortices ( Fromang & Nelso"nll2005h and zonal 
flows dJohans en et al. 2009). 

We also observe that the energy cascade is dominated by the 
magnetic and exchange fluxes down to the resistive scale, the ki- 
netic flux and shear fluxes being almost negligible. Note also that 
the shear fluxes are always positive. This is due to the anisotropy 
of MRI turbulence in which shearing waves have statistically a 
larger amplitude than leading waves (see section [T2b . 

In the Pm = 0.0625 case, the kinetic flux becomes domi- 
nant at subresistive scale {K > 10), indicating that the cascade 
becomes essentially hydrodynamic below the resistive scale, 
as expected for small Pm MHD turbulence (see also Fig. |6]l. 
Moreover, the kinetic flux dominates the kinetic shear transfer 
term at least for the larger k, which indicates that, as far as 
the non linear transfers are concerned, the cascade is close to 
isotropic at small scales, as noted in section 13.21 Finally, note 
that none of the fluxes reach a plateau at intermediate scales, 
which would be expected in the presence of an inertial range. 
This indicates that the k^^^- kinetic spectrum found in Fig. [1] is 
not properly speaking an inertial spectrum. 

3.5. Energy transfer locality 

To test the locality of the energy transfer in spectral space, we 
have plotted the transfer functions in the Pm - 0.25 case for 
several values of K: at the injection scale (K - 1), in the in- 
termediate range (K = 5) and in the resistive range (K = 20). 
We first plot the kinetic to kinetic and magnetic to magnetic 
transfers on Fig. I?] The transfers Tw and Tt,/, obtained at all 
scales perfectly illustrate local energy exchanges. Energy is 
taken from wavenumbers slightly smaller than K and is trans- 
ferred to wavenumbers slightly larger than K, except (not sur- 
prisingly) for K = 1 . As expected from the energy flux, we also 
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Fig. 6. Zoom on the energy fluxes in the dissipative range for 
Pm - 0.0625 (log-log representation with K > 10). 

find that the cascade is direct, with energy going from small to 
large wavenumbers. Finally, note that the T,,,, transfers are al- 
ways much smaller than the Tth transfers above the dissipation 
range, illustrating the magnetically dominated energy transfer 
described above. 

We next plot the exchange transfers Tyi, and Thy on Fig. [8] 
We note in this case that the scales involved in each transfer 
are much broader In particular, Tyb measured at the resistive 
scale (K = 20) has contribution coming from all scales, includ- 
ing the largest injection scales. This effect can also be seen in 
Tbv(Q,K = 1), which exhibit a very long tail toward large k, 
down to the resistive scale. Comparing directly these transfers to 
Tyy and Tbb show that these terms are highly non local. In turns, 
this indicates that the exchange flux computed in the previous 
section is non local. The results in the Pm - 0.0625 case are not 
shown here as they are very similar to the Pm = 0.25 case. 

We note that despite the non locality of the energy ex- 
changes, the overall cascade direction is still forward, confirm- 
ing our previous interpretation regarding the exchange flux. We 
also remark that the shear transfer terms are local by definition 
as they transfer energy to neighbouring shells. 

4. Discussion 

In this paper, we have described some properties of the turbu- 
lent cascade found in incompressible MRI turbulence. We have 
shown that compared to isotropic MHD turbulence, the presence 
of a mean shear led to several new transfer terms and introduced 
a source of anisotropy. We have computed the effect of each non 
linear term and found that all the terms contribute to a direct 
cascade of energy (from large to small scales) but some terms 
involved non local transfers in Fourier space. This non locality 
is due to the Lorentz force and to the magnetic stretching term of 
the induction equation (combined here in the exchange transfer 
term). We have also shown that magnetic helicity, although non 
zero, was totally negligible and should not play any role in the 
behaviour of MRI turbulence. 

The presence of non local transfer terms in the MRI turbu- 
lent cascade is the most important finding of this work. It in- 
dicates that in principle, the large scales - responsible for the 
transport - can directly interact with the small dissipative ones 
through non linear terms. This direct interaction could of course 
explain the correlation observed between Pm and the turbulent 
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Fig. 7. Transfert function T^XQ, K) and Thb(Q, K) in the Pm - 0.25 run for K - 1 ; 5; 20. These transfers are local in Fourier space 
(see text). 




transport of angular momentum a (iLongaretti & Lesurll2010h . 
However, it should be pointed out that nonlocal transfe rs were 
alrea dy found in isotropic MHD turbulence (e.g. Alexakis et alj 
12007 ). Therefore, MRI turbulence has nothing special regarding 
the nature of these non linear transfers. 

Although some nonlinear terms are found to directly connect 
injection and dissipation scales in current simulations, one might 
wonder if this would be true in a more realistic setup where the 
injection and dissipation scales are separated by a wide range 
of scales (typically 10'"). In other words, what is the maximum 
scale separation thes e terms can connect? A partial answer to this 
question is given bv lAluie & EvinkI (j2010). In order to describe 
their result, let us define the structure functions: 

6vi^p = {\v(x + 1)- v{x)\P) (24) 

In the inertial range the structure function depends only on |Z| and 
6vip oc /fp, where (p is the structure function index of order p. It 
is then possible to derive an upper bound to the nonlinear transfer 
terms thanks to the Holder inequality. Applying this procedure 
to the non local transfer Tub, Aluie & Evink (.2010.) found 

\Tub{Q, K)\ < iconst)Q'-^""K-^^'" (25) 

where Q and K are dyadic (octave) wavenumbers and K > Q/2. 
Similar terms can be obtained for r^,, and K < Q/2. If one as- 
sumes Iroshnikov-Kraichnan theory, one has - - 3/4. On 
the contrary, considering Goldreich-Sridhar (GS) phenomenol- 
ogy, which should be valid for MRI turbulence, one gets ^^'n = 

^3II =3/2 and = = 1. In all these cases, ( |25] | indicates 
that the non-locality of these transfer terms cannot extend over 
several decades, with a typical scaling Tub{Q,K) ~ e{K/Qy^^^ 
for GS turbulence (e being the usual turbulence energy injection 
rate). 



We therefore conclude that the non locality in Fourier space 
is somewhat relative. Although Tub and Tbu are non local com- 
pared to Tuu and Tbb, these terms should be local when one 
considers transfers over several decades. Unfortunately, sepa- 
rating the injection scale from the dissipative scales by several 
decades is numerically difficult. It is even harder for MRI tur- 
bulence since the injection term is rather broad in spectral space 
compared to forced turbulence. Assuming the injection and dis- 
sipation scales both spread over one decade in Fourier space, 
one typically needs 20, 000^ simulations to get a 2 decades iner- 
tial range in which non local transfers are significantly reduced. 
This kind of resolution is for the moment out of reach of the best 
computational facilities. 

Nevertheless, we can conjecture that if the Pm-a correlation 
is effectively due to the non local transfers, then it should van- 
ish when the injection and dissipation scales are well separated 
as it is the case in some accretion discs. Although this conclu- 
sion looks rather reassuring for the relevance of today simula- 
tions regarding small scale dissipation, it does not tell us what 
the asymptotic value of a is in this limit, nor how MRI turbu- 
lence behaves when the scale separation is not achieved, a situa- 
tion which probably occurs in the inner regions of protoplanetary 
discs where A,j is not very large. 



Appendix A: Shearing waves approach to the shear 
transfer term 

The shell-filter decomposition can be properly defined using a 
projector operator II on the field F in the sheared frame: 

[n^,(F)](x',f)= 2 Fk'(t)exp{ik' -x') (A.l) 
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Fig. A.l. Evolution of shearing waves in the presence of a fixed 
shell in Fourier space. Some waves can either enter or exit the 
shell as time evolves (see text). 



where I,j is the shell containing all the shearing waves with a 
norm between Kj - SK/2 and Kj + 6K/2: 

= i^Kj - 6K/2 < \k' + ^Qfe;,feJ < Kj + 6K/2^ (A.2) 

Our notation indicates that the projected F is function of space 
and time. Note also that [njf^(F)] is real for real fields F. 

As it can be seen, the waves included in 2y(f) change with 
time. This is to be expected as shearing waves should in princi- 
ple enter and exit shells as the one defined above (see Fig. lA.lb . 
As a result, the projector operator rij^^ has an explicit time de- 
pendence which leads to non trivial transfer effects. The above 
projector operator can be written using Heaviside function 0: 

[UK^(F)](x',t) = Yj®{\k(t)\-Kj + 6K/2^ 

k' 

X 0^ - \k(t)\ + Kj + 6K/2^ 



Fk'{t)exp(ik' ■ x'^ . 



(A3) 



where we have defined k(t) as a function of k' as in ( fTTT l — ( fTlT l. 
One next defines the energy within a shell = (D^ (F))/2 

where {•) denotes a volume average. The energy equation then 
reads 



^ 

dt 



= {nKj(F) 



dUK^jF) 
dt 



) 



= (n^,(F)%(5,F)) + (n^,(F)[5,%](F)) (A.4) 

The first term on the right handside leads to the terms obtained 
in isotropic turbulence (Alexakis et al. 2007) and introduced in 
Eqs. ( fTsT i and ( fTSI l along with the injection terms ly K and Ih^K- 
The second one however is due to shearing waves entering and 
exiting the shells. Using (IA.3b . it is possible to obtain an exact 
(though singular) expression of the operator time derivative: 



5,[%](F) 



qQ.kykx{t) 



6i^k{t)\ - Kj + 5KI2)j -si^- \k{f)\ + Kj + 6KI2)j 



Fk'(t)exp(ik' ■ x'^ . 



(A.5) 



This expression can be interpreted easily. As an example, let us 
consider waves with kyk^(t) > 0. Then, the first 6 function repre- 
sents waves entering the shell, the second delta represents waves 
leaving the same shell, and the factor in front of the delta func- 
tions quantifies the "flux" of waves going through a shell bound- 
ary. This interpretation is similar to the phenomenological pic- 
ture one can have of waves traveling through a fixed shell in the 
unsheared Fourier space (Fig. lA.lt . 
We then deduce from (IA.4b 



dt 



<%(F)%(5,F)) 

q£lkyk_,(t) Fl,Fk' 



(A.6) 



m\ 



x\6{\k{t)\ - Kj + 6K/2^ - 5^ - \k(t)\ + Kj + 6K/2^ 

where we have used the propertjQ 6(x)6(x) - 6(x)/2. As ex- 
pected, this expression shows two contributions to the energy 
evolution inside a shell: a volume contribution and a surface 
contribution equal to the energy of the waves entering and leav- 
ing the shell. Introducing the equation of motions in sheared 
space (|2l) — (O in the above relation leads to the energy equa- 
tions ([T5]l — ( fT6b . 

Note finally that using the relation 



d,\k(t)\ 



qQ.kykAt) 



m\ 

one can write (IA.6b as 



dE^ 
dt 



F* F ' 



{A.l) 



(A.8) 



where tk' is the instant when the wave k' enters or exits the shell 
Kj and ew - +1 for an entering/exiting wave. This somewhat 
simpler expression has the same interpretation as ( IA.6b . 

Appendix B: Unsheared Fourier transform 
approach to the shear transfer term 

It is possible to understand the origin of the shear transfer term 
(IA.5b starting from the equations of motion in unsheared coor- 
dinates ([T]i — (|2]i with an appropriate use of continuous Fourier 
transforms. Shear-periodic functions are not absolutely inte- 
grable, but this difficulty can be circumvented because their con- 
tinuous Fourier transform is well-defined as a distribution. To 
demonstrate this point, let us consider a 2D infinite medium in 
which a field F obeys the model equation 



d,F(x, t) - qQxdyF{x, t) = 0. 

Let us introduce the unsheared Fourier transform: 



(B.l) 



F(k, t) 
F(x 



(2^ III ■ 
, f) - JJJ dk F{k, t) exp(/fe ■ a;). 



(B.2) 



\k{t)\ 



This can be shown differentiating the relation 0^(x) = @(x) 
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If F obeys the shearing sheet boundary conditions, the solu- 
tion of the equation is a Fourier series Fj, of the form ( fTOb with 
Fourier coefficients Ff,/ . Its Fourier transform in the unsheared 
spectral space is then: 



Fp(k,t) = y6(k-kQ(t))Fk' 



(B.3) 



where fco(0 = k' + Vk't and where we have defined the 
"Fourier velocity" 14 - q^K^x, evaluated in k' . By construc- 
tion Fp{k, t) is solution of: 

d,F{k) + qQ.kydkF{k) = 0, (B.4) 

which is the Fourier transform of our model equatiorQ. The 
left handside of this equation can be interpreted as a comoving 
derivative of F{k, t) in Fourier space with the Fourier velocity Vk ■ 
This equation tells us that the amplitude of the waves is constant 
when one moves at velocity Vk in Fourier space consistently with 
the form of the solution Fp. It can also be interpreted as constant 
amplitude shearing waves, as expected. 

It is then possible to introduce the projector operator, now 
time-independent as it is defined in unsheared coordinates: 

[nK,{F)]{x, t) = ^ dk @{\k\ - Kj + 6K/2j 
X - + Kj + 

X F(k,t)exp(ik-x). (B.5) 
{Ul (F)}/2) time varia- 



As previously, the shell energy {E^ 
tion follows from: 



dE 



= (n^,(F)5,[n^,(F)]) 

= -(%(F)%(Vfc ■ (VkF))) 



(B.6) 



where ( IB.4l i has been used in the last equality as well as ■ 
Vk - 0. Because the volume average selects the zero frequency 
contributions, this leads us back to ( IA.6I 1 with the help of the 
relation 



dk,&{\k\-C)dkF(k) 



6(\k\ - C)F(k). 



and with the use of 6(x)6(x) - 6(x)/2 and (IB. 31 ). 

In the spirit of the integral expressions used in this appendix, 
Eq. (IA.6b can be recast in integral form by introducing the energy 
density in Fourier space &'^{k, t) defined by 



&\k,t)^Yj^{k-ko{t))—. 



Fk'Fl, 



k' 

With this definition, 
E^, = JJJ dk &{\k\ - K-)&(- 

while 

dE'' 



I -I- K^) S''(k, f) 



= JJJ dkVk■n6''{k,t)[6{\k\-K-)-6i- 



dkVk■n &^{k, t) 



(B.7) 
(B.8) 

(B.9) 
(B.IO) 



In these relations, — K + 6K/2 has been defined; the 
second integration is performed on the surface of the shell Kj = 
dKj and n is the normal to this surface. The last expression has 
an explicit flux form. 

This approach can be applied to the original shearing sheet 
MHD equations; the time dependence due to the shear term will 
produce the shear flux contribution just computed. 
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